TIL contexture predicts checkpoint blockade response in melanoma patients
In this case study, we will use ProjecTILs to interpret human scRNA-seq T cell data in the context of a murine TIL atlas. We are going to use the single-cell dataset generated by Sade-Feldman GSE120575 to illustrate interspecies mapping of human cells on a stable murine atlas using gene orthologs.
Background
It remains unclear why some patients respond to checkpoint blockade therapy while others do not.
In the study by Sade-Feldman et al. (2018) Cell, the authors characterize transcriptional profiles of immune cells from melanoma patients before and after immune checkpoint blockade, with the goal to identify factors that associate with success or failure of checkpoint therapy. They reported that the balance between two CD8 T cell states found in tumor tissue (essentially, memory-like vs exhausted-like) is linked to tumor regression following checkpoint blockade. And in particular, that the frequency of TCF7+ CD8+ T cells in tumors predicts response and better survival.
The single-cell expression data GSE120575 consists of CD45+ single cells from 48 tumor biopsies of melanoma patients before or after checkpoint inhibitors treatment, and sequenced using the Smart-seq2 protocol. Meta-data on response to therapy (Responder vs. Non-responder) is also available on the same GEO identifier.
Here we use the development version of ProjecTILs >0.4.1, which implements ortholog gene conversion between human and mouse.
R Environment
Check & load R packages
Sys.setenv(R_REMOTES_NO_ERRORS_FROM_WARNINGS = "true")
if (!requireNamespace("renv")) install.packages("renv")
library(renv)
renv::restore()
# Load development version of ProjecTILs implementig human-mouse ortholog
# conversion
remotes::install_git("https://gitlab.unil.ch/carmona/ProjecTILs.git", ref = "v0.4.1")
library(Seurat)
library(ProjecTILs)
library(gridExtra)scRNA-seq data preparation
Download the count matrix and metadata from Gene Expression Omnibus (GEO), and store as Seurat object.
cached.object <- "SadeFeldman.seurat.rds"
if (!file.exists(cached.object)) {
library(GEOquery)
geo_acc <- "GSE120575"
datadir <- "input/SadeFeldman"
gse <- getGEO(geo_acc)
series <- paste0(geo_acc, "_series_matrix.txt.gz")
system(paste0("mkdir -p ", datadir))
getGEOSuppFiles(geo_acc, baseDir = datadir)
## Load expression matrix and metadata
exp.mat <- read.delim(sprintf("%s/%s/GSE120575_Sade_Feldman_melanoma_single_cells_TPM_GEO.txt.gz",
datadir, geo_acc), header = F, sep = "\t")
genes <- exp.mat[c(-1, -2), 1]
cells <- as.vector(t(exp.mat[1, 2:16292]))
samples <- as.factor(t(exp.mat[2, 2:16292]))
exp.mat <- exp.mat[c(-1, -2), 2:16292]
colnames(exp.mat) <- cells
rownames(exp.mat) <- genes
meta <- read.delim(sprintf("%s/%s/GSE120575_patient_ID_single_cells.txt.gz",
datadir, geo_acc), header = T, sep = "\t", skip = 19, nrows = 16291)
meta <- meta[, 1:7]
treat <- factor(ifelse(grepl("Post", samples), "Post", "Pre"))
response <- factor(meta$characteristics..response)
therapy <- factor(meta$characteristics..therapy)
## Create Seurat object and add meta data
query.object <- CreateSeuratObject(counts = exp.mat, project = "SadeFeldman",
min.cells = 10, min.features = 50)
rm(exp.mat)
query.object@meta.data$Sample <- samples
query.object@meta.data$Time <- treat
query.object@meta.data$Response <- response
query.object@meta.data$Therapy <- therapy
saveRDS(query.object, file = cached.object)
} else {
query.object <- readRDS(cached.object)
}Some basic statistics - cells per group (Pre vs. Post, Therapy, Responder vs. Non-responder)
Post Pre
10363 5928
anti-CTLA4 anti-CTLA4+PD1 anti-PD1
517 4121 11653
Non-responder Responder
10727 5564
Post_P1 Post_P1_2 Post_P10
291 353 308
Post_P10_T_enriched Post_P11 Post_P12
88 313 323
Post_P13 Post_P13_T_enriched Post_P14
265 184 329
Post_P15 Post_P16 Post_P16_T_enriched
304 369 94
Post_P17 Post_P17_myeloid_enriched Post_P18
238 43 357
Post_P19 Post_P19_myeloid_enriched Post_P2
328 75 301
Post_P2_T_enriched Post_P20 Post_P21
79 358 213
Post_P22 Post_P23 Post_P23_2
260 358 362
Post_P28 Post_P28_2 Post_P3
369 372 359
Post_P3_2 Post_P30 Post_P4
358 369 421
Post_P4_T_enriched Post_P5 Post_P5_2
254 292 336
Post_P6 Post_P6_T_enriched Post_P7
308 92 231
Post_P8 Post_P8_T_enriched Pre_P1
327 82 229
Pre_P12 Pre_P15 Pre_P2
330 304 337
Pre_P20 Pre_P24 Pre_P25
323 338 371
Pre_P26 Pre_P27 Pre_P28
333 343 337
Pre_P29 Pre_P3 Pre_P31
163 245 351
Pre_P33 Pre_P35 Pre_P4
452 238 311
Pre_P6 Pre_P7 Pre_P8
288 368 267 Select only baseline (Pre-treatment) samples
Post_P1 Post_P1_2 Post_P10
0 0 0
Post_P10_T_enriched Post_P11 Post_P12
0 0 0
Post_P13 Post_P13_T_enriched Post_P14
0 0 0
Post_P15 Post_P16 Post_P16_T_enriched
0 0 0
Post_P17 Post_P17_myeloid_enriched Post_P18
0 0 0
Post_P19 Post_P19_myeloid_enriched Post_P2
0 0 0
Post_P2_T_enriched Post_P20 Post_P21
0 0 0
Post_P22 Post_P23 Post_P23_2
0 0 0
Post_P28 Post_P28_2 Post_P3
0 0 0
Post_P3_2 Post_P30 Post_P4
0 0 0
Post_P4_T_enriched Post_P5 Post_P5_2
0 0 0
Post_P6 Post_P6_T_enriched Post_P7
0 0 0
Post_P8 Post_P8_T_enriched Pre_P1
0 0 229
Pre_P12 Pre_P15 Pre_P2
330 304 337
Pre_P20 Pre_P24 Pre_P25
323 338 371
Pre_P26 Pre_P27 Pre_P28
333 343 337
Pre_P29 Pre_P3 Pre_P31
163 245 351
Pre_P33 Pre_P35 Pre_P4
452 238 311
Pre_P6 Pre_P7 Pre_P8
288 368 267 ProjecTILs
Load reference TIL atlas - if it’s not present in the working directory, it will be downloaded from the repository
[1] "Loading Default Reference Atlas..."
[1] "/Users/admin/gitHub/ProjecTILs_CaseStudies/ref_TILAtlas_mouse_v1.rds"
[1] "Loaded Reference map ref_TILAtlas_mouse_v1"Run ProjecTILs projection algorithm - version 0.4.1 supports human ortholog conversion (human.ortho=T)
[1] "Using assay RNA for query"
[1] "1685 out of 5928 ( 28 % ) non-pure T cells removed. Use filter.cells=FALSE to avoid pre-filtering (NOT RECOMMENDED)"
[1] "Transforming expression matrix into space of mouse orthologs"
[1] "Aligning query to reference map for batch-correction..."
Projecting corrected query onto Reference PCA space
[1] "Projecting corrected query onto Reference UMAP space"Plot global projection of human TIL data over the reference in UMAP space (non-T cell transcriptomes have been automatically removed)
Predict the cell states in the query set
query.projected <- cellstate.predict(ref=ref, query=query.projected)
table(query.projected$functional.cluster)
CD4_NaiveLike CD8_EarlyActiv CD8_EffectorMemory CD8_NaiveLike
82 413 789 1071
CD8_Tex CD8_Tpex Tfh Th1
583 188 170 612
Treg
335 Interestingly, expression of marker genes of projected human TILs correspond fairly well to those of the murine reference atlas: e.g. Pdcd1, Havcr2 and Entpd1 expression in terminally exhausted CD8 TILs (CD8_Tex), co-expression of Pdcd1, Tox, and Tcf7 (at a modest level) in the CD8_Tpex state with a higher expression of Ifng; high expression of Tcf7, Ccr7 with low expression of Pdcd1 or cytotoxic molecules in the Naive-like states (that might include central memory cells); Cxcr5 and Tox coexpression in follicular helper T cells; Foxp3, Tox, Havcr2, Entpd1 in regulatory CD4 T cells, etc.
query.list <- SplitObject(query.projected, split.by = "Response")
plot.states.radar(ref, query = query.list, min.cells = 50, genes4radar = c("Foxp3",
"Cd4", "Cd8a", "Tcf7", "Ccr7", "Gzmb", "Pdcd1", "Havcr2", "Tox", "Entpd1", "Cxcr5",
"Ifng", "Cxcl13", "Xcl1", "Itgae"))
Note that projections and comparisons are performed in the ortholog space of murine genes - to check the names of human-mouse orthologs you can examine the conversion table for genes of interest:
data(Hs2Mm.convert.table)
which.genes <- c("TCF7","GZMB","CD8B","PDCD1","ITGAE")
Hs2Mm.convert.table[Hs2Mm.convert.table$Gene.HS %in% which.genes, ] Gene.stable.ID.HS Gene.HS Gene.MM Alt.symbol
6099 ENSG00000083457 ITGAE Itgae Itgae
7381 ENSG00000254126 CD8B Cd8b1 Cd8b1
7835 ENSG00000188389 PDCD1 Pdcd1 Pdcd1
11096 ENSG00000081059 TCF7 Tcf7 Tcf7
16688 ENSG00000100453 GZMB Gzmb Gzmd,Gzme,Gzmf,Gzmc,Gzmn,Gzmg,Gzmb
Alt.symbol.HS
6099 ITGAE
7381 CD8B2,CD8B
7835 PDCD1
11096 TCF7
16688 GZMBResponse to therapy
Now to the interesting part. Let’s visualize the projection and TIL contexture of tumors that responded or not responded following immune checkpoint blockade:
query.list <- SplitObject(query.projected, split.by = "Response")
pll <- list()
pll[[1]] <- plot.projection(ref, query.list[["Responder"]]) + ggtitle("Responder")
pll[[2]] <- plot.statepred.composition(ref, query.list[["Responder"]], metric = "Percent") +
ggtitle("Responder") + ylim(0, 40)
pll[[3]] <- plot.projection(ref, query.list[["Non-responder"]]) + ggtitle("Non-responder")
pll[[4]] <- plot.statepred.composition(ref, query.list[["Non-responder"]], metric = "Percent") +
ggtitle("Non-responder") + ylim(0, 40)
grid.arrange(grobs = pll, ncol = 2, nrow = 2, widths = c(1.5, 1))
In non-responders, there is a clear enrichment in the terminally exhausted CD8 state (CD8_Tex), while responders are enriched in Naive-like states.
To better examine differences in TIL contexture between responders and non-responders, we can visualize the fold-change of T cell state frequency between the two groups:
which.types <- table(query.projected$functional.cluster) > 20
stateColors_func <- c("#edbe2a", "#A58AFF", "#53B400", "#F8766D", "#00B6EB", "#d1cfcc",
"#FF0000", "#87f6a5", "#e812dd")
states_all <- levels(factor(ref$functional.cluster))
names(stateColors_func) <- states_all
cols_use <- stateColors_func[names(which.types)][which.types]
# Responder vs non Responder
query.list <- SplitObject(query.projected, split.by = "Response")
norm.c <- table(query.list[["Non-responder"]]$functional.cluster)/sum(table(query.list[["Non-responder"]]$functional.cluster))
norm.q <- table(query.list[["Responder"]]$functional.cluster)/sum(table(query.list[["Responder"]]$functional.cluster))
foldchange <- norm.q[which.types]/norm.c[which.types]
foldchange <- sort(foldchange, decreasing = T)
tb.m <- melt(foldchange)
colnames(tb.m) <- c("Cell_state", "Fold_change")
pll <- list()
ggplot(tb.m, aes(x = Cell_state, y = Fold_change, fill = Cell_state)) + geom_bar(stat = "identity") +
scale_fill_manual(values = cols_use) + geom_hline(yintercept = 1) + scale_y_continuous(trans = "log2") +
theme(axis.text.x = element_blank(), legend.position = "left") + ggtitle("Responder vs. Non-responder")
Indeed, CD4 and CD8 Naive-like states are the most enriched in responders compared to non-responders, while terminally exhausted Entpd1+ Havcr2+ CD8 TILs are the most under-represented, confirming the observation of the original paper that a higher frequency of TCF7+ CD8 TILs is associated with response to immune checkpoint therapy.
Conclusions
Taking advantage of the ortholog mapping functionality of ProjecTILs, we have illustrated how to effortlessly analyze human scRNA-seq data in the context of a reference murine TIL atlas. Gene expression profiles confirmed that T cells are accurately projected in major CD4+ and CD8+ categories, as well as in more specific subtypes (CD8_Tex, CD8_Tpex, naive-like, Follicular helper, Th1, and regulatory CD4 cells).
Comparison of transcriptional profiles and cell states of TILs at baseline from responding vs non-responding melanoma patients, confirmed the original observation by Sade-Feldman et al that the frequency of TCF7+ CD8 TILs correlates with checkpoint therapy responsiveness.
However, the TCF7+ CD8 TIL population associated to response seems to correspond to a naive-like state and not to the (PD-1+ TOX+) Precursor exhausted state characterized in murine cancer and chronic infection models (Siddiqui et al. 2019; Miller et al. 2019).
Moreover, this naive-like TIL population is unlikely to be tumor-specific, especially in the absence of strong evidence for clonal expansion, as the frequency of tumor-reactive cells among PD-1- CD8 TILs has been shown to be very low in melanoma tumors (Gros et al. 2014).
This might seem at odds with other studies showing that the presence of PD-1+ ‘exhausted-like’ CD8 T cells are predictive of response to checkpoint blockade in melanoma and non-small cell lung cancer (Daud et al. 2016; Thommen et al. 2018). However, in this dataset most tumors from both responders and non-responders did contain a pool of exhausted-like TILs expressing PD-1 (PDCD1), CD39 (ENTPD1), CXCL13, CD103 (ITGAE), even if its frequency was higher among non-responders.
Therefore, the presence of a pool of PD-1+ CD8 T cells in the tumor might be required, but not sufficient, for therapeutic success. Other factors, including total amount of CD8 TIL infiltration, spatial distribution, etc. might be correlated with the presence of naive-like CD8 TILs and with improved response. For example, we also observed that compared to non-responding, responding tumors had higher frequencies of Th1-like and naive-like CD4 T cells and lower of regulatory CD4 T cells, which might also contribute to the reduced anti-tumor response following immunotherapy.
ProjecTILs analysis in the context of a stable atlas gives a more complete picture of the T cell states that are associated with immunotherapy response.
Further reading
Dataset original publication - Sade-Feldman et al. (2018) Cell
ProjecTILs case studies - INDEX - Repository
References
Daud, A.I., Loo, K., Pauli, M.L., Sanchez-Rodriguez, R., Sandoval, P.M., Taravati, K., Tsai, K., Nosrati, A., Nardo, L., Alvarado, M.D., Algazi, A.P., Pampaloni, M.H., Lobach, I.V., Hwang, J., Pierce, R.H., Gratz, I.K., Krummel, M.F., Rosenblum, M.D., 2016. Tumor immune profiling predicts response to anti-PD-1 therapy in human melanoma. J. Clin. Invest. 126, 3447–3452. https://doi.org/10.1172/JCI87324
Gros, A., Robbins, P.F., Yao, X., Li, Y.F., Turcotte, S., Tran, E., Wunderlich, J.R., Mixon, A., Farid, S., Dudley, M.E., Hanada, K., Almeida, J.R., Darko, S., Douek, D.C., Yang, J.C., Rosenberg, S. a, 2014. PD-1 identifies the patient-specific in filtrating human tumors. J. Clin. Invest. 124, 2246–59. https://doi.org/10.1172/JCI73639.2246
Miller, B.C., Sen, D.R., Al Abosy, R., Bi, K., Virkud, Y.V., LaFleur, M.W., Yates, K.B., Lako, A., Felt, K., Naik, G.S., Manos, M., Gjini, E., Kuchroo, J.R., Ishizuka, J.J., Collier, J.L., Griffin, G.K., Maleri, S., Comstock, D.E., Weiss, S.A., Brown, F.D., Panda, A., Zimmer, M.D., Manguso, R.T., Hodi, F.S., Rodig, S.J., Sharpe, A.H., Haining, W.N., 2019. Subsets of exhausted CD8+ T cells differentially mediate tumor control and respond to checkpoint blockade. Nat. Immunol. 20, 326–336. https://doi.org/10.1038/s41590-019-0312-6
Siddiqui, I., Schaeuble, K., Chennupati, V., Fuertes Marraco, S.A., Calderon-Copete, S., Pais Ferreira, D., Carmona, S.J., Scarpellino, L., Gfeller, D., Pradervand, S., Luther, S.A., Speiser, D.E., Held, W., 2019. Intratumoral Tcf1+PD-1+CD8+ T Cells with Stem-like Properties Promote Tumor Control in Response to Vaccination and Checkpoint Blockade Immunotherapy. Immunity 50, 195–211.e10. https://doi.org/10.1016/J.IMMUNI.2018.12.021
Thommen, D.S., Koelzer, V.H., Herzig, P., Roller, A., Trefny, M., Dimeloe, S., Kiialainen, A., Hanhart, J., Schill, C., Hess, C., Prince, S.S., Wiese, M., Lardinois, D., Ho, P.C., Klein, C., Karanikas, V., Mertz, K.D., Schumacher, T.N., Zippelius, A., 2018. A transcriptionally and functionally distinct pd-1 + cd8 + t cell pool with predictive potential in non-small-cell lung cancer treated with pd-1 blockade. Nat. Med. 24, 994. https://doi.org/10.1038/s41591-018-0057-z